function query_z = gridInterp2d(raster, x0, y0, dx, dy, query_x, query_y)
    query_x = (query_x - x0) / dx + 1;
    query_y = (query_y - y0) / dy + 1;
    [size_y, size_x] = size(raster);
    [xx, yy] = meshgrid(1:size_x, 1:size_y);
    query_z = interp2(xx, yy, raster, query_x, query_y, "cubic", nan);
end
